Learning how to classify raster data

# load libraries
library(raster)
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.1-10, (SVN revision 622)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-3

Import CHM (Canopy Height Model)

# import TEAK CHM
chm <- raster("../NEONdata/D17-California/TEAK/2013/lidar/TEAK_lidarCHM.tif")

# test plot
plot(chm, main = "Plot TEAK CHM using raster package")

# test stretch image
image(log(chm), main = "Plot TEAK CHM image \nJust render pixels that stretch to fill the space")

Deal with 0 values

# plot hist to identify 0 values
hist(chm)
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 32% of the raster cells were used. 100000 values used.

# set 0 values to NA
chm[chm == 0] <- NA

# look at hist again
hist(chm, xlab = "Tree height (m)")  # lots of tall trees

Read in aspect layer

# read in aspect
aspect <- raster("../NEONdata/D17-California/TEAK/2013/lidar/TEAK_lidarAspect.tif")

# plot aspect
plot(aspect, main = "Plot of TEAK aspect data")

Let’s look at how aspect relates to tree height

Are trees taller on north or south aspects?

First, we will set north-facing slopes to 1 and south-facing slopes to 2. * North-facing slopes are 0-45 and 315-360 degrees, Class = 1 * South-facing slopes are 135-225 degrees, Class = 2

Create classification matrix

# create array of values
class.m <- c(0, 45, 1,
             45, 135, NA,
             135, 225, 2,
             225, 315, NA,
             315, 360, 1)

# shape values into matrix
rcl.m <- matrix(class.m,
                ncol = 3,
                byrow = TRUE)

# let's look at our matrix
rcl.m
##      [,1] [,2] [,3]
## [1,]    0   45    1
## [2,]   45  135   NA
## [3,]  135  225    2
## [4,]  225  315   NA
## [5,]  315  360    1

Now reclassify the raster

# reclassify raster using matrix
asp.ns <- reclassify(aspect,
                     rcl.m)

# plot reclassified raster
plot(asp.ns, main = "North- and south-facing slopes")

Export Geotiff

Save reclassified matrix as output.

writeRaster(asp.ns,
            file = "../outputs/TEAK/TEAK_nsAspect2.tif",
            options = "COMPRESS=LZW", 
            NAflag = -9999)

# compress to keep file size as small as possible

Create a mask

We want to create a raster mask that will compare two rasters only at the selected points of north- and south-facing slopes. This masks out the data we are not interested in.

Important: When using a mask to compare two different rasters, need to pay attention to: Are the two rasters the same resolution? Are the locations matched up or is there georeference data in the raster?

# create a mask from asp.ns
asp.ns
## class       : RasterLayer 
## dimensions  : 577, 543, 313311  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 325963, 326506, 4102905, 4103482  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 0, 2  (min, max)
# import NDVI layer
ndvi <- raster("../NEONdata/D17-California/TEAK/2013/spectrometer/veg_index/TEAK_NDVI.tif")

# plot it
plot(ndvi,
     main = "NDVI for TEAK site")

# mask data
nsFacing.ndvi <- mask(ndvi, asp.ns)

# plot it
plot(nsFacing.ndvi, 
     main = "NDVI for north- and south-facing slopes \nTEAK field site")

Subset aspect raster to only select north-facing slopes

# subset raster so only north-facing slopes have a value
asp.n <- asp.ns == 1
asp.n[asp.n == 0] <- NA

# test plot
plot(asp.n,
     main = "North-facing slopes \nTEAK field site")

Use this new north-facing raster to mask NDVI

# mask NDVI with north-facing slopes
nFacing.ndvi <- mask(ndvi, asp.n)

# plot
plot(nFacing.ndvi,
     main = "NDVI for north-facing slopes \nTEAK field site")

Subset aspect raster to only select south-facing slopes

# subset raster so only south-facing slopes have a value
asp.s <- asp.ns == 2
asp.s[asp.s == 0] <- NA

# test plot
plot(asp.s,
     main = "South-facing slopes \nTEAK field site")

Use this new south-facing raster to mask NDVI

# mask NDVI with south-facing slopes
sFacing.ndvi <- mask(ndvi, asp.s)

# plot
plot(sFacing.ndvi,
     main = "NDVI for south-facing slopes \nTEAK field site")